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was  conducted  to  correlate  calculational  approaches  with  data. 
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CONVERSION  TABLE 


Conversion  factors  for  U.S.  customary  to  metric  (SI)  units  of  measurement 


To  Convert  From 

To 

Multiply 

angstrom 

meters  (m) 

1.000  OOOXE-10 

atmosphere  (normal) 

kilo  pascal  (kPa) 

1.013  25XE+2 

bar 

kilo  pascal  (kPa) 

1.000  OOOXE+2 

barn 

metcr^  (m^) 

1.000  OOOXE-28 

British  Thermal  unit  (thcrmochemlcal) 

joule  (J) 

1.054  350XE+3 

calorie  (thermochemical) 

joule  (J) 

4.184  000 

cal  (thermochemlcal)/cm^ 

mega  joule /m^(MJ/m*) 

4.184  OOOXE-2 

curie 

glga  becquerel  (GBq)* 

3.700  OOOXE+1 

degree  (angle) 

radian  (rad) 

1.745  329XE-2 

degree  Fahrenheit 

degree  kelvln  (K) 

tK=(t®f  + 459.67)/ 1.8 

electron  volt 

joule  (J) 

1.602  19X£>-19 

erg 

joule  (J) 

1.000  000  X  E-7 

erg/ second 

watt  (W) 

1.000  000  X  E-7 

foot 

meter  (m) 

3.048  OOOXE-1 

foot-pound-force 

joule  (J) 

1.355  818 

gallon  (U.S.  liquid) 

meter^  (m^) 

3.785  412  XE-3 

inch 

meter  (m) 

2.540  000  X  E-2 

Jerk 

joule  (J) 

1.000  000  X  E+9 

joule /kilogram  (J/Kg)  (radiation  dose 
absorbed) 

Gray  (Gy) 

1.000  000 

kilotons 

terajoules 

4.183 

kip  (1000  IbO 

newton  (N) 

4.448  222  X  E+3 

klp/lnch^  (ksl) 

kilo  pascal  (kPa) 

6.894  757  X  E+3 

ktap 

newton-second /m^  (N-s/m^) 

1.000  000  X  E+2 

micron 

meter  (m) 

1.000  000  X  E-6 

mil 

meter  (m) 

2.540  000  X  E-5 

mile  (International) 

meter  (m) 

1.609  344  X  E+3 

ounce 

kilogram  (kg) 

2.834  952  X  E-2 

pound-force  (Ibf  avoirdupois) 

newton  (N) 

4.448  222 

pound-force  inch 

newton-meter  (N-m) 

1.129  848X  E-1 

pound-force  /  Inch 

newton /meter  (N/m) 

1.751  268  X  E+2 

pound-force  /  foot^ 

kilo  pascal  (kPa) 

4.788  026  X  E-2 

pound-force /inch^  (psi) 

kilo  pascal  (kPa) 

6.894  757 

pound-mass  (Ibm  avoirdupois) 

kilogram  (kg) 

4.535  924  X  E-1 

pound-mass-foot^  (moment  of  inertia) 

kilogram-meter^  (kg-m^) 

4.214  011  X  E-2 

pound-mass  /  foot^ 

kilogram /meter^  (kg/m^) 

1.601  846XE+1 

rad  (radiation  dose  absorbed) 

Gray  (Gy)^* 

1.000  000  X  E-2 

roentgen 

coulomb /kilogram  (C/kg) 

2.579  760  X  E-4 

shake 

second  (s) 

1.000  OOOX  E-8 

slug 

kilogram  (kg) 

1.459  390X  E+1 

torr  (mm  Hg,  O^C) 

kilo  pascal  (kPa) 

1.333  22  X  E-1 

•The  becquerel  (Bq)  is  the  SI  unit  of  radioactivity;  Bp  =  1  event/s. 
••The  Gray  (Gy)  Is  the  SI  unit  of  absorbed  radiation. 
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SECTION  1 


INTRODUCTION 

1.1  BACKGROUND  AND  OBJECTIVE. 

This  technical  report  describes  computational  modeling  of  underground 
tunnel  response  in  support  of  the  Underground  Technology  Program  (UTP)  which 
is  a  multi-year  investigation  Into  the  vulnerability  of  underground  structures.  The 
overall  program  includes  computational  modeling,  material  modeling,  laboratory 
testing,  and  field  testing  to  improve  the  ability  to  predict  the  response  and  failure 
of  underground  structures  subjected  to  ground  shock  due  to  near-surface 
explosions.  The  emphasis  is  on  deeply  buried  tunnels  with  little  or  no 
reinforcement. 


1.2  SCOPE  AND  SUMMARY. 

The  computational  modeling  effort  involved  three  major  areas  of 
investigation:  1)  benchmark  activity,  2)  parametric  study  and  3)  laboratory  test 
simulation.  The  benchmark  activity  was  a  step-by-step  series  of  idealized 
problems  which  addressed  various  aspects  of  tunnel  response  in  jointed  rock 
masses.  The  parametric  study  considered  the  systematic  sensitivity  of  loading 
environment,  material  characterization  and  geometric  conditions  as  well  as 
computational  approaches  on  tunnel  response.  Numerical  simulation  of  the  SRI 
International  lab-scale  HE  experiments  on  tunnel  response  in  intact  limestone 
was  conducted  to  correlate  calculational  approaches  with  data. 

As  one  of  several  participants  in  the  computational  modeling  effort,  it  is 
clear  that  a  step-by-step  procedure  of  material  model  verification,  free-field 
loading  environment  definition  and  parametric  response  calculations  need  to  be 
performed  in  order  to  reduce  uncertainties  and  improve  the  ability  to  predict  the 
response  and  failure  of  underground  structures  subjected  to  ground  shock  due  to 
near-surface  explosions. 
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SECTION  2 


BENCHMARK  ACTIVITY 

Various  numerical  approaches  to  problems  of  tunnel  dynamics  are 
compared  with  each  other  and,  wherever  possible,  with  exact  analytic  solutions 
by  Logicon  RDA  [Simons,  1992].  The  medium  is  an  idealization  of  a  jointed  rock 
mass.  The  intact  rock  is  linear  elastic-plastic  with  a  pressure  dependent  failure 
surface  and  associated  plastic  flow  law.  There  are  two  orthogonal  sets  of  equally 
spaced  joints.  Each  joint  is  nonlinear  elastic  In  the  normal  direction  and  linear 
elastic  with  Coulomb  friction  in  shear. 

The  TRT  approach  utilizes  the  EXCALIBUR  finite  element  method  [Ito, 
England  &  Nelson,  1981]  to  represent  the  intact  rock  and  joints  with  two  different 
types  of  models  for  jointed  rock,  an  explicit  one  where  the  joints  are  treated 
separately  (in  the  near  field  of  the  tunnel),  and  an  implicit  one  where  their 
properties  are  lumped  together  with  those  of  the  intact  rock  (in  the  far  field).  The 
TRT  implicit  joint  model  admits  arbitrary  constitutive  behavior  in  both  the  intact 
rock  and  joint,  and  by  enforcing  internal  compatibility  and  stress  equilibrium 
derives  a  super-element  representing  the  combined  deformation  due  to  both 
joints  and  intact  rock. 


2.1  PROBLEMS  WITH  ANALYTICAL  SOLUTIONS. 

The  first  five  problems  are  quasi-static  driven  by  boundary  displacements 
which  are  consistent  with  homogeneous  (uniform)  strain  throughout  the  region  of 
interest.  This  does  not  mean  that  the  actual  strains  will  be  uniform;  if  there  are 
joints  then  the  actual  strains  will  not  be  uniform.  But  in  fact  the  stresses  and 
strains  in  the  intact  rock  and  joint  material  will  separately  be  homogeneous  at 
each  point  of  the  imposed  strain  path.  This  opens  the  possibility  of  direct 
analytical  solution  of  these  problems  for  comparison  with  numerical  results. 

Another  way  of  viewing  the  situation  is  that  each  of  these  problems 
reduces  to  nothing  more  than  finding  the  response  of  a  single  implicit  element 
around  a  specified  strain  path.  This  Is  strictly  a  material  response  question. 
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equations  of  motion  or  compatibilty  among  elements  play  no  role  whatsoever.  It 
is  interesting  to  note  that  TRT  was  the  only  participant  to  produce  single-element 
solutions  to  all  of  the  first  five  problems. 

2.2  TWO-DIMENSIONAL  PROBLEMS. 

Simplicity  makes  anaiytic  solutions  possible.  In  contrast  the  two- 
dimensional  problems  have  fields  which  vary  both  spatially  and  temporally  so 
that  analytic  solutions  are  not  generally  feasible.  The  free-fleld  problem  concerns 
deformations  of  a  wedge-shaped  section  of  an  annulus  in  plane  strain,  as  shown 
In  Rgure  2-1,  with  the  entire  region  containing  vertically  and  horizontally  jointed 
rock.  The  top  edge  (inner  arc)  Is  loaded  with  the  pressure  pulse  shown  in  the 
figure,  while  shear  tractions  are  zero.  The  left  and  right  sides  have  roller 
boundaries,  making  the  left  side  a  plane  of  symmetry  (the  right  side  is  not, 
because  the  effective  anisotropy  due  to  jointing  makes  the  material  unsymmetric 
about  that  plane).  The  lower  edge  (the  outer  arc)  is  a  transmitting  boundary. 

Most  of  the  region  Is  to  be  modeled  implicitly,  except  for  a  rectangular 
region  extending  2.5  tunnel  diameters  (12.5  m)  in  all  directions  from  the  on-axis 
point  at  R  =  500  m.  This  region  Is  modeled  explicitly,  in  anticipation  of  the  tunnel 
situated  there  in  the  final  problem.  The  free-field  problem  is  essential  to  clearly 
define  the  loading  environment  at  the  tunnel  range.  The  geometry  and  loading  in 
the  final  problem  are  shown  in  Figure  2-2,  which  are  precisely  the  same  as  the 
free-field  problem  but  now  there  is  a  lined  tunnel  in  the  center  of  the  rock  island. 

Based  on  a  combination  of  physical  understanding  of  wave  propagation 
and  material  behavior  and  comparison  with  the  analytic  solutions,  three  of  the 
numerical  approaches,  including  TRT,  are  judged  [Simons,  1992]  to  have 
produced  credible  results  to  the  final  problem  of  a  lined  tunnel  in  jointed  rock 
mass  engulfed  by  a  cylindrically  divergent  stress  wave. 
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SECTION  3 


PARAMETRIC  STUDY 


The  parameter  study  provides  a  systematic  investigation  of  the  sensitivity 
of  tunnel  response  to  a  variety  of  factors; 

-  Geometric  Considerations 

-  range  from  source 

-  diameter 

-  divergence 

-Loading  Environment 

-  peak  incident  stress 

-  pulse  duration 

-  Material  Characterization 

-  dilatancy 

-  extension/compression  strength  ratio 

-  flow  rule 

-  Numerical  Scheme 

-code 

g 

The  scope  of  the  study  (initially)  involves  nine  independent  parameters  or  2  = 
512  two-dimensional  calculations  using  a  wedge-shaped  computational  model 
and  linear  elastic-plastic  material  model  similar  to  that  in  the  benchmark  activity 
but  without  joints.  The  participants  are  Weidlinger  Associates  (WA)  and  TRT. 
The  calculational  data  from  the  parameter  study  is  analyzed  by  Logicon  RDA 
[Pucik  &  Curry,  1993]  using  design-of-experiments  methodologies  [Box,  Hunter  & 
Hunter,  1978]. 
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3.1  PROBLEM  DEFINITION. 


The  initial  problem  definition  for  geometry,  material  properties,  loading  and 
calculation  matrix  are  given  In  Tables  3-1  to  3-4,  respectively. 


3.2  CONCLUSIONS. 

The  key  conclusions  from  this  parameter  study  are  as  follows: 

-  Tunnel  Failure  Sensitive  to  Three  Main  Parameters 

-  peak  incident  stress 

-  extension/compression  strength  ratio 

-  flow  rule 

-  Crush  Strength  Not  A  Very  Sensitive  Parameter 

-  Numerical  Approach  (Calculator)  An  Important  Parameter 

-  Geometric  Scaling  Applies 

-  range  and  diameter  combined  into  single  parameter 

-  if  and  only  if  numerical  grid  scaled  geometrically 
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Table  3-1.  Geometry. 


Table  3-2.  Material  Properties. 


Table  3-3.  Loading. 


Parameter 

Value(s)  I 

Peak  Incident  stress 

Positive  phase  duration  of 
velocity  relative  to  DUG-1  C 

f/f+DUG 

1 

0.5 

Wavefront  divergence 

0 

1  1 

Table  3-4.  Calculation  Matrix. 


Run 

R 

(m) 

(MPa) 

(deg) 

m 

1 

500 

200 

25 

1 

10 

1 

133.3 

2 

500 

200 

25 

0.5 

5 

0.5 

133.3 

3 

500 

200 

0 

1 

5 

0.5 

66.7 

4 

500 

200 

0 

0.5 

10 

1 

66.7 

5 

500 

100 

25 

1 

10 

0.5 

66.7 

6 

500 

100 

25 

0.5 

5 

1 

66.7 

7 

500 

100 

0 

1 

5 

1 

133.3 

8 

500 

100 

0 

0.5 

10 

0.5 

133.3 

9 

250 

200 

25 

1 

5 

1 

66.7 

10 

250 

200 

25 

0.5 

10 

0.5 

66.7 

11 

250 

200 

0 

1 

10 

0.5 

133.3 

12 

250 

200 

0 

0.5 

5 

1 

133.3 

13 

250 

100 

25 

1 

5 

0.5 

133.3 

14 

250 

100 

25 

0.5 

10 

1 

133.3 

15 

250 

100 

0 

1 

10 

1 

66.7 

16 

250 

100 

0 

0.5 

5 

0.5 

66.7 

SECTION  4 


LABORATORY  TEST  SIMULATIONS 

SRI  International  performed  both  static  tunnel  tests  [Simons  et  al,  1 993] 
and  spherical  wave  tunnel  (SWAT)  tests  [Klopp  et  al,  1993]  in  iimestone  blocks  to 
compare  the  measured  loading  environment  and  tunnel  closures  with  those  from 
various  computational  models.  The  TRT  approach  utilizes  relatively  fine  zoning 
to  capture  the  free-field  environment  with  ten  circumferential  elements  across  the 
tunnel  radius  and  radial  elements  one-fifth  of  the  radius  near  the  tunnel,  while  the 
liner  is  a  simple  membrane.  Both  two-invariant  and  three-invariant  failure  models 
are  developed  for  the  limestone  based  on  the  RE/SPEC  material  properties 
[Possum,  1993].  All  numerical  simulations  were  performed  ‘pre-tesf  before  the 
measured  data  were  made  available.  In  addition,  preliminary  calculations  were 
conducted  to  verify  the  computational  approach  and  material  modeling. 

4.1  STATIC  TUNNEL  TESTS. 

Three  static  tunnel  tests  In  limestone  are  simulated,  called  ST1,  ST2,  and 
ST3.  The  limestone  specimens  are  cylinders  30.5  cm  (12  In.)  in  diameter  with  a 
tunnel  located  at  the  specimen  midheight.  The  height  of  the  specimen  for  ST1  Is 
30.5  cm  (12  in.)  and  the  tunnel  Is  unlined  and  has  a  diameter  of  19.1  mm  (0.750 
in.).  For  ST2  and  ST3  the  specimen  height  is  45.8  cm  (18  in.)  and  the  tunnels 
are  lined  (fully  annealed  3003  aluminum  tube  with  70  MPa  strength).  The 
specimens  are  loaded  by  two  Independent  hydraulic  pressures,  a  vertical  load 
applied  to  the  top  and  bottom  surfaces  of  the  cylinder  and  a  confining  load 
around  the  outer  surface  of  the  cylinder.  For  ST  1  and  ST3  the  loading  history  is 
a  single  cycle  of  loading  (191  MPa  and  139  MPa,  respectively)  and  unloading. 
For  ST2  the  specimen  is  loaded  to  125  MPa,  unloaded,  reloaded  to  156  MPa, 
and  then  unloaded. 

The  closure  comparisons  of  the  numerical  simulations  of  ST1  using  two- 
invariant  (M-S)  and  three-invariant  (W-W)  strength  models  based  on  RE/SPEC- 
93  properties  show  little  difference  at  the  crown-invert  (Figure  4-1)  but  a  factor  of 
three  increase  at  the  springline  (Figure  4-2)  due  to  strength  reduction  (70%)  in 
TXE  ,  which  is  more  consistent  with  the  measured  closure  data.  Numerical 
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simulations  of  ST2  produce  similar  results  as  ST1  since  the  confining  conditions 
are  essentially  uniaxial  strain  with  no  shear  failure  in  the  far-field  at  the  loading 
boundary.  In  contract,  both  the  numerical  simulations  and  laboratory 
experiments  of  STS  are  very  sensitive  and  produce  a  wide  variation  in  tunnnel 
closure  due  to  the  divergent  loading  which  can  induce  shear  failure  in  the  whole 
specimen. 

4.2  SPHERICAL  WAVE  TUNNEL  TESTS. 

The  spherical  wave  tunnel  (SWAT)  tests  involve  a  charge  of  PETN 
detonated  in  a  block  of  limestone  containing  tunnels.  The  numerical  simulations 
involve  the  two  initial/boundary  value  problems  given  in  Figure  4-3.  The 
problems  are  identical  except  for  the  location  of  the  tunnel,  either  14.5  cm  or  19.5 
cm  from  the  center  of  the  charge.  The  measured  particle  velocity  history  shown 
in  the  figure  is  applied  to  the  inside  of  the  7.50-cm  radius  cavity. 

A  rate-enhanced  three-invariant  (W-W)  failure  with  symmetric  crush  (Cap) 
model  consistent  with  RE/SPEC  and  ARA  static  (10'®/s)  data,  ARA  extrapolated 
dynamic  (10’^^  -  10‘"^/s)  strength  data  and  LLNL  shock  hugoniot  (10'^/s)  crush 
data  is  developed  for  numerical  simulation  of  SWAT.  Assuming  log-linear 
interpolation  with  no  SWAT  free-field  data,  this  uncalibrated  rate-enhanced 
model  gives  about  a  factor  of  two  increase  in  crush  strength  (initial  cap  location) 
and  33%  increase  in  shear  strength. 

The  rate-enhanced  model  is  also  calibrated  using  log-bilinear  interpolation 
to  match  SWAT  (10"*^^^  -  10'^®/s)  free-field  environment  at  both  the  particle 
velocity  stations  (Figure  4-4)  as  well  as  the  stress  at  the  range  of  the  tunnels 
(Figure  4-5).  The  comparison  of  the  calibrated  rate-enhanced  response  to  the 
measured  data  at  the  14.5-cm  tunnel  is  8%  versus  7%  crown-invert  closure 
(Figure  4-6)  and  4%  versus  2%  springline  closure  (Figure  4-7).  This  is  a  very 
good  comparison  given  the  uncertainty  in  strain-rate  effects  on  shear  strength 
including  extension/compression  strength  ratio.  It  should  be  noted  that  while  the 
rate  independent  model  appears  to  give  better  correlation  with  measured 
springline  closure  (Figure  4-7),  this  is  very  misleading  since  this  model  does  not 
reproduce  the  measured  free-field  loading  environment  (Figures  4-4  and  4-5). 
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Figure  4-2.  Springline  closure  -  ST  1 . 
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SECTION  5 


RECOMMENDATIONS 

As  one  of  several  participants  in  the  computational  modeling  effort, 
it  is  clear  that  a  step-by-step  procedure  of  material  model  verification,  free-field 
loading  environment  definition  and  parametric  response  calculations  need  to  be 
performed  to  reduce  uncertainties  and  improve  the  ability  to  predict  the  response 
and  failure  of  underground  structures  subjected  to  ground  shock  due  to  near¬ 
surface  explosions.  A  major  uncertainty  Is  strain-rate  effects  on  shear  strength 
including  extension/compression  strength  ratio  of  rocks  of  interest. 
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